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OO I Brownian Dynamics algorithms are widely used for simulating soft-matter and biochemical 
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systems. In recent times, their application has been extended to the simulation of coarse- 
grained models of cellular networks in simple organisms. In these models, components move 
by diffusion, and can react with one another upon contact. However, when reactions are 
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^ I incorporated into a Brownian Dynamics algorithm, attention must be paid to avoid violations 

I of the detailed-balance rule, and therefore introducing systematic errors in the simulation. 

We present a Brownian Dynamics algorithm for reaction-diffusion systems that rigorously 
obeys detailed balance for equilibrium reactions. By comparing the simulation results to 
exact analytical results for a bimolecular reaction, we show that the algorithm correctly 
reproduces both equilibrium and dynamical quantities. We apply our scheme to a "push- 
pull" network in which two antagonistic enzymes covalently modify a substrate. Our results 
highlight that the diffusive behaviour of the reacting species can reduce the gain of the 
' response curve of this network 

o ■ 
cr: 

I. INTRODUCTION 

' Most, if not all, biological processes are regulated by biomolecules, such as proteins and DNA, 



which chemically and physically interact with one another in what are called biochemical networks. 
These networks are often highly non-linear, which means that mathematical modelling is critical 
I for understanding and predicting their behaviour. The dominant paradigm has been to consider 

' the living cell to be a spatially homogeneous environment, analogous to a well-stirred reactor. It 

. is increasingly recognised, however, that the cell is a highly inhomogeneous environment, in which 

I compartmentalisation, scaffolding and localised interactions are actively exploited to enhance the 

regulatory function of biochemical networks. This means that it becomes important to not only 
describe the biochemical network in time, but also in space. 

In this manuscript, we present an algorithm for simulating biochemical networks in time and 
' space that is based upon Brownian Dynamics. Brownian Dynamics is a stochastic dynamics scheme, 

in which the solvent is treated implicitly; only solutes are described explicitly. The forces expe- 
rienced by the solutes contain a contribution from the interactions with the other solutes and 
a random part, which is the dynamical remnant of the collisions with the solvent molecules. A 
pioneering BD algorithm was introduced by Ermak and McCammon [l|, and detailed, atomistic 
Brownian Dynamics simulations have beenperformed to study the dynamics of enzyme-substrate 
and protein-protein association reactions [2 [3, 0, [B, 0, Q, Si- More recently, Brownian Dynamics 
has not only been used to stud y th e association between two proteins, but also to simulate networks 
of interacting biomolecules 

In order to simulate these large systems at the biologically 
relevant length and time scales, molecules are coarse-grained to the level of simple geometrical 
objects, which can diffuse and react with other chemical species in a confined geometry. 

While Brownian Dynamics algorithms for simulating biochemical networks are based upon a 
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simplified description of the molecules and their interactions, they do go beyond the conventional 



kinetic Monte Carlo schemes to simulate biochemical networks [12f |. These algorithms are based 
upon the zero-dimensional chemical master equation, and, as such, they take into account the 
discrete nature of the components and the stochastic character of their interactions. However, 
they assume that at each instant in time the particles are uniformly distributed in space. In 
contrast, Brownian Dynamics based algorithms take into account not only the particulate nature 
of the molecules and the probabilistic character of their interactions, but also that at any moment in 
time the particles may be non-uniformly distributed in space. Brownian Dynamics thus accounts 
for both temporal and spatial fluctuations of the components. Moreover, it allows for spatial 
gradients and for localised interactions in the network. 

Recently, a number of stochastic techniques have been developed that make it possible to 
simulate biochemical networks in time and space. Some techniques are based upon Brownian 
Dynamics while others are based upon the reaction-diffusion master equation 0, Q 

[lil . [l^ . The advantage of Brownian Dynamics based techniques is that they are truly particle- 
based, which means that they do not have to rely on a mesoscopic length and time scale on 
which the system is well-stirred. We have recently developed an entirely novel algorithm, called 
Green's Function Reaction Dynamics (GFRD), which is particle-based scheme to simulate biological 
networks in time and space, like Brownian Dynamics. However, in contrast to Brownian Dynamics, 
GFRD is an event-driven algorithm, which uses Green's functions to concatenate the propagation 
of the particles in space with the chemical reactions between them. This makes GFRD orders of 
magnitude more efficient than Brownian Dynamics when the concentrations are below 0.1 — 1/iM. 
For higher concentrations, or for reactions near surfaces, brute-force Brownian Dynamics is more 
efficient, because of the smaller computational overhead per time step. 

Although the main idea of applying Brownian Dynamics to reaction-diffusion systems is straight- 
forward, a number of ingredients has to be examined with care. One is: to which processes do the 
association and dissociation rates as used in the simulations correspond to? To the intrinsic rates, 
which are the reaction rates at contact, or to the effective rates that also take into account the effect 
of diffusion? The other issue is detailed balance. Biochemical networks often contain reactions that 
do not consume energy or are otherwise driven out of equilibrium. These equilibrium reactions 
should obey detailed balance. Even though a number of Brownian Dynamics based algorithms have 
been presented 13, U], this question has, to our knowledge, not been systematically addressed. 

In this paper, we present a Brownian Dynamics algorithm that rigorously obeys detailed balance 
and is thus able to reproduce the equilibrium properties of a reaction-diffusion system. In Section iHl 
we derive our algorithm on the basis of the statistical mechanics of chemical reactions. The 
algorithm is subjected to stringent tests in Section Hill besides equilibrium properties, we test also 
how well the algorithm reproduces the dynamical behavior of a bimolecular reaction, for different 
values of the time step At. A comparison with a stochastic algorithm that does not account for 
spatial fluctuations is also presented. Finally, in Section IIVI we show an illustrative application 
of our algorithm to a simple coarse-grained model of a chemical species under the action of two 
enzymes operating in opposite directions (the so-called "push-pull" model system). Simulations 
conducted with our BD algorithm show that both spatial and temporal fluctuations reduce the 
gain of the response of the system. 
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II. METHODS 
A. Detailed balance 

Before we present the outline of the algorithm in the next section, we discuss the detailed- 
balance rule that must be obeyed for equilibrium reactions. To this end, we will consider the 
elementary bimolecular reaction: 

A + B^C {kon,kos). (1) 

Here A;on is the macroscopic forward rate for the association of molecules A and B, and /coff is the 
macroscopic backward rate for their dissociation. The macroscopic expression for the equilibrium 
constant for this reaction is 

^ _ kon_ _ [C] , . 

~ [A][By 

where [X] is the concentration of the species X. 



In a spatially-resolved model, we can decompose the reaction ([T]) in two steps |l7l |: 



A + B A-B^C. (3) 

^D.b fed 

In the first step, particles A and B find each other and form an encounter complex A ■ B, which 
has not yet reacted to a final product; this occurs via a diffusion-limited rate kj) = AttRD, where 
R = {Ra + Rb)/'^ is the cross section with Rx the diameter of particle X, and D = Da + Db, 



with Dx is the diffusion constant of species X 13] . Given that the particles are in contact, the 



reaction can then proceed according to the intrinsic reaction rate fca. The rates and A:D,b denote 
the intrinsic dissociation rate and the rate at which the particles in the encounter complex diffuse 



away into the bulk, respectively [17|, [iSl]. It can be shown [171] that the equilibrium constant is 
given by 

~k,- feoff 

and that the macroscopic forward and backward rate constants are given by, respectively, 

— - — + — (5) 

^on 

koS fed ' 

We will use Brownian Dynamics to simulate not only the diffusive motion of the particles, but 
also the reactions between them. At each step of the algorithm, each particle is given a trial 
displacement according to a distribution that follows from the diffusion equation, as described 
below. If the move does not lead to an overlap with another particle, the move is accepted. 
Importantly, this procedure naturally simulates the formation of the encounter complex with a rate 
feD, provided that the step sizes are smaller than the diameters of the particles. If two particles 
are close to each other, and thus form an encounter complex, a trial displacement of one of the 
two can lead to an overlap; this overlap leads to a reaction with a probability, as derived below, 
that is consistent with the intrinsic reaction rate fea. Conversely, at each step of the algorithm a 
product particle C can dissociate with a probability consistent with the intrinsic dissociation rate 
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k^. If a trial dissociation move is accepted, then the particles A and B have to be put back in the 
encounter complex. The question, however, is: at which distance should the particles be put back 
relative to each other? The Brownian Dynamics scheme makes an error in the dynamics of order 
At. This might suggest that the precise location is not critically important, as long as the distance 
is smaller than Ar ~ \/DAt. However, not every choice obeys detailed balance, and, as we will 
show, a choice that does not obey detailed balance will lead to systematic errors. We now derive 
the detailed balance condition. 

The detailed-balance condition for one given pair of particles A and B is 

-^unbound (l")'^ri'u_»l3 — -Pbound-fb^U) (7) 

where -Pbound is the probability that the two particles are bound, and -Punbound(r)t^r is the 
probability that the particles A and B are separated by a vector between r and r + dr. We 
now first derive the ratio -Pbound/(-funbound(r)dr). To this end, let us consider the probability 
P{r^'',r^'',r^,^;{NA,NB,Nc})dr^^dr^Bdr^c that the system has {Na,Nb,Nc) molecules and 

that these molecules are located at positions {r^, • • • ,r^''^}, {r]^, • • • ,r^^}, {r^, • • • ,r^'^}. This 
probability is given by 

P{r'^\4^,r^^;{NA,NB,Nc}) = PN{NA,NB,Nc)xr{4^4^,r^'=^\{NA,NB,Nc}), (8) 

where Pi\f{NA, Nb, Nc) is the probability that the system has {Na, Nb, Nc) molecules and V is 
the conditional probability density that a given number {Na, Nb, Nc} of molecules occupy those 
particular positions. As discussed in more detail in Appendix|Al P]\[{Na, Nb, Nq) is given by 

Na Nb Nc yNA+Ng+Nc -. 

Pn{Na,Nb,Nc) = NaINbINcI Q' 

where Q is the partition function of the system. The conditional probability density is the proba- 
bility density of finding {Na, Nb, Nc} indistinguishable ideal particles in a volume V: 

Vi4^r-^, 4o\{NA,NB,Na}) = (10) 

Combining the above equations, yields the following expression for the probability density: 

Na Nb Nc 

T-,f Na N„ Nr tat at at t\ ^A,cmlB,cni1c ,cm /-I -I \ 

Pi^A^'^B^'^C '^{^A,Nb,Nc}) = Q , (11) 

where qx,cm is the partition function corresponding to the degrees of freedom od the center of 
mass of a particle X, and Q is canonical partition function of the system. The ratio between the 
probability densities of being in a state after and before the transition is: 



P{4^-\r^^-\r^.^'+';{NA-l,NB-l,Nc + l}) _ qc,cm 



P{r^^ , r^^ , r^^ ; {Na, Nb, Nc}) qA,cmqB,cm 

By taking Na = 1, Nb = l,Nc = 0, we obtain Pbound/(-Punbound(r)c^r): 

Abound P{rc;0,0,l)drc qC, cm 

^'unbound(r)c^r P(rA, r^; 1, 1, 0)dr^(irs qA,cmqB,cradr dr ' 

Using Eq. HI the detailed-balance condition, Eq. \7\ then becomes 

-fbound ^u^b 

^'unbound(r)cir Ph^udr k^dr 



(12) 



(13) 



(14) 
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As discussed above, the association between the particles in the encounter comlpex to form the 
product C consists of a two-step process: 1) a "generation" move or "trial" move, in which an 
overlap is generated with a probability -Pgcn,f) 2) an "acceptance" move, in which the overlap is 
accepted with probability Pacc.f! the product of the probabilities of these moves is related to the 
intrinsic reaction rate ka_. Similarly, the dissociation move also consists of two steps: 1) a "trial" 
move, in which the dissociated particles are put at a vector between r and r + dr with probability 
-fgen,b(r)dr; 2) an "acceptance" move, in which the trial move is accepted with probability -Paccb! 
the product of the probabilities of these moves is related to the intrinsic dissociation rate constant 
/cd- The detailed-balance condition can thus be written as [l^ 

-fbound -fgcn.f (r)-Pacc,f f^K'\ 

-Punbound(r)dr -Pgcn,bdr(r)Pacc,b ^d^r ' 

This is the principal result of this section. Below, we discuss in detail how this rule is implemented 
in our BD scheme. In appendix [Aj we discuss how this detailed-balance rule is related to the 
detailed-balance rule for a well-stirred system, where we do not account for the positions of the 
particles in space. 



B. Simulation scheme 



It is instructive to consider the association between one particle A and one particle B. We can 
assume without loss of generality that = 0, i.e. that the A particle does not diffuse in the 
simulation box. It is then convenient to position it at the center of the box. The single B particle 
moves by free diffusion with diffusion coefficient Db = D. At every simulation step, the system is 
propagated by a fixed time At. 

In the absence of the A particle, the motion of the B particle is simply described by the Einstein 
equation: 

d 

— p{r',t + At\r,t) = DV^p{r',t + At\r,t), (16) 

where p{r' , t + At\r, t) is the probability of finding the particle at position r' at time t + At, given 
that it was at r at time t. We know with certainty the position of the particle at the initial time. 
We also impose that at time t + At the probability of finding the particle in space vanishes as we 
move far away from the initial position r. We can then formulate the following boundary conditions 
for Eq. ([16]): 

p(r',t + At|r,t) = 5(r'-r), (17) 



p(|r'| ^ oo,t + At|r,t) = 0. (18) 

The solution of (|16p with conditions (jl7p and (jlSp is a Gaussian function, whose variance is pro- 
portional to At: 

1 ( (r' - r)2 1 

p(r', t + Atlr, t) = -—7:7 exp <^ - ^ —j- } . (19) 

' ' ' ^ (2-2DAt)3/2 ^\ 2 -20 At] ^ ' 

This time-dependent probability distribution can be used to generate new positions for the B 
particle at every time step At [20|. 

In the presence of the A particle, a reaction can occur when the B particle overlaps with the A 
particle. In order to describe the association and dissociation reactions, we have to specify Pgen(r), 
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-facc,fi -fgen,bi ^acc,b in such a Way that detaield balance, Eq. [151 is obeyed. We first discuss 
Pgen,f(r), then the two quantities related to the backward move, Pgen,b(r) and -Pacc,b) and then we 
discuss the probability by which the trial association move (the overlap) should be accepted, -Pacc,f ■ 
The quantity Pgen,f(r) can be computed analytically: let us consider the single particle A 
held fixed in a center of a large box, whose edges lie far enough to be neglected in the following 
derivation. Using a polar reference frame whose origin coincides with the center of the A sphere, we 
can compute the probability that a B particle initially at position r is displaced to a position r' G E, 
where E is the excluded volume for B (a sphere, centered in the origin, with radius R = Ra+Rb)'- 

p(r^S)=/ r'^dr smOdO I dipp{r' ,t + At\r,t) = g{r, At). (20) 
Jo Jo Jo 

The function g{r) can be computed analytically, is radially symmetric and depends on the Brownian 
Dynamics time step At. Details are given in Appendix [Bj We will not indicate anymore the 
dependence of various quantities on At, since this parameter is kept constant during the whole 
simulation. We remind the reader that in the function g{r), r represents the position from which the 
B particle leaves, given that the move led to an overlap with A. We set then Pgen, f (r) =5(r)r2(0, (/?), 
where VL[9, ip) is the uniform angular distribution on the sphere. 

Dissociation is modelled as a first order reaction event, with a Poissonian distribution of waiting 
times: P{t) = kdexp(—kdt). The probability that the reaction has not happened at time t is 
then S{t) = 1 — Jq P{t')dt' = exp(— /cdt). Therefore, the probability a reaction does happen is 
1 — exp(— fcdi) — k^At if k^At ^ 1. If we choose time steps At such that At ^ 1/A;d, the 
probability that an event happens within At can then be approximated to k([At. We therefore 
accept the dissociation move with a probability Pace, b = k^^At. 

Once we have determined that a dissociation event has happened, we must determine a new 
position for the B particle in the reaction box. The crux of our BD algorithm is to generate a 
reverse move according to a probability distribution Pgen,b(r) that is the same as that by which the 
forward move is generated, Pgen,f(r) = g{r)fl{6,4>), but properly renormalised. The normalisation 
factor can be obtained by integrating Pgen,f(r) over all initial distances r: 

/"OO f 

J dr J dn n{e, ip)g{r, At)r'^ = 47rI(At) < V, (21) 

where 1 = g{r)r'^dr. During a dissociation move, the particle is thus put at a vector between r 
and r + dr according to Pgcn,h{^)dr = ■^g{r)Q,{9, (f). 

Using Eq. [T5l we can now obtain the desired acceptance probability for the forward move: 

p -f bound -fgen,b p 

^unbound gen,f 

_ A;a g{r)n{9,ip)dr ^ 



fcddr g{r)0,{9,Lp) AttI 
k^At 



(22) 



The above expression has a meaningful interpretation: the intrinsic association rate k^. can be 
written as the product of two factors: 1) a collision frequency 47r//At and 2) the probability 
-facc,f that a collision leads to a reaction. The dominant contributions to the integral / come 
from distances r that are short compared to V DAt (see Eqs. [20l and [2T]) . In the limit that the 
time step At 0, the rate /ca should thus approach the intrinsic association rate fca as used in 
theories of diffusion- influenced reactions [TtI]; here the intrinsic association rate fca is defined as 
the association rate given that the particles A and B are in contact. We also note that this is the 



intrinsic association rate as used in Green's Function Reaction Dynamics 2l|, [2^ 
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C. Algorithm outline 

Let us consider a system with M particles of type B and one particle of type A, held fixed at 
the center of box of volume V. For convenience, we choose as initial state the situation in which 
there is no bound state C. 

1. Generate an initial position for the B particles in the available volume. 

2. Select randomly one of the particles among species B and C. 

3. (a) If the particle is type B, for each cartesian coordinate, generate a new position according 

to a Gaussian distribution with zero mecin and stctiidard deviation y/^iD^t'. 3^new — 
a;old+^^(0, V2DAi), where At the Br ownian Dynamics time step. 

(b) If the displacement move leads to an overlap of the B particle with A, that is if |r^ — 
r^l < Ra + Rb, attempt a reaction according to a probability -Pacc,f = feaAt/(47rI). 

(c) If the trial reaction move is accepted, remove the B particle from the box, and substitute 
the A particle with a C. This new particle is not diffusing in the box. 

(d) If the trial reaction move is rejected, put the B particle back to its original position. 

4. (a) If the particle is type C, try a backward reaction with probability Pace, b = kdAt. 

(b) If the trial reaction move is accepted, substitute the C particle with an A particle, 
create a new B particle whose radial position is drawn from the normalised distribution 
g{r)/I and the angular position from the uniform distribution Q,{9, (p). If this leads to 
an overlap with another B particle, reject the move. 

(c) If the trial reaction move is rejected, keep the identities and positions of the particles. 

5. Repeat step 2. and 3. or 4. M times, then increase the simulation time by At. 

Keeping particle A and C fixed could mimic for example a system where one reactant is anchored 
to some rigid scaffold. A relevant biological example is the binding of proteins to DNA in a bacterial 
cell, particle A representing a binding site on the DNA, typically in proximity of some gene. In 
this case, the motion of A is only related to the fluctuations of the polymer, which happen on time 
scales much longer than the diff^usion of proteins in the bacterial cytoplasm, and can therefore be 
neglected. The scheme could straightforwardly be extended to the situation in which the A particle 
also moves, or cases with more reaction channels. 



III. TESTS 

In this Section, we check the BD scheme by comparing the simulation results with analytical 
results. Our scheme was built upon a series of assumptions, which should all be satisfied simulta- 
neously. In particular, 1) the time steps should not be too large, as a BD algorithm is not able 
to resolve the system at time scales below the time step At; 2) the acceptance probabilities for 
the forward and backward reactions should be small (-Pacc,f ^ 1) -facc,b ^ !)• In particular, the 
algorithm does not resolve the precise moment in time when the association and dissociation events 
happen. Dynamical quantities could therefore exhibit systematic errors, which must vanish in the 
limit At^O. On the other hand, on long enough time scales, even dynamical poperties should be 
reproduced, provided that the conditions listed above are obeyed. In the case of two particles, the 
probability distribution p(r,t|ro) of finding the particles separated by a vector r at time t given 
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that initially they were separated by fq, as well as their survival probability 5(t|ro) [17|, have been 
computed analytically [23]; they will provide a stringent test for the dynamics of our scheme. As 
our algorithm obeys detailed balance, equilibrium quantities, such as the average time spent in the 
bound state, must be correctly reproduced for all time steps At. 



A. Irreversible Reactions 



We begin by simulating the irreversible reaction A + B — ^ C, within the following setup: a 
single particle A is held fixed in an unbounded system, and a single particle B is positioned on a 
spherical surface at an initial distance ro from A, with a random angle. The particles have the same 
radius: Ra = Rb = R/2. We run the algorithm for a time tsim, and we record the final radial position 
of the particle B. In the case that a reactive event happens before tsinn we stop the run. After a 
large number of runs, we collect the final positions of the B particle in an histogram, normalised 
to the fraction of B particles which have survived until the final time. This histogram should 
reproduce the irreversible probability distribution Pirr{r,tsim\ro,0) [S]- This quantity represents 
the probability of finding the two particles at time tsim separated by a distance r, given an initial 
separation of tq at to = 0- We note that this probability distribution is not normalised to unity: 
the integral over space of pi^r is the survival probability of the particle, which is the probability 
that the particle has not reacted at the final time. Formally: 

/•oo 

4vr / pi„{r,t\ro)r'^ dr = Sirr{t\ro). (23) 
JR 

We are thus able to simultaneously test our algorithm twice: comparing the analytical curve with 
the profile of our histogram, and the area of the histogram with the analytical value of the survival 
probability. 

Results are collected in Figure [D we simulate the irreversible reaction for 4 different simulation 
times, from tsim = 10~^r to tsim = 10~^t, where t = R?/D is the natural time scale of the system. 
Particles are initially positioned at contact: tq = R. We set the time step At = 10~^tsim) which 
corresponds to Pacc.f < 0.14. It is seen that both the shape and the area of the irreversible 
probability disitribution function is correctly captured by our algorithm. In the case of tsim = 10~^r, 
however, we needed to use At = 10~^tsim- In the Inset of Fig. [H we show that in this last case, 
larger time steps lead our BD scheme to underestimate the survival probability. The deviation 
from the analytical results for large At is due to the interplay bewteen a number of assumptions. 
One is that 1 — exp(Pacc,f) — -facc,f- A more important factor is that we compare the numerical 
results against analytical results of an analysis in which kon corresponds to the intrinsic association 
rate for two particles that are at contact [23], while in our scheme the particles can already react 
when they are separated by a distance ~ \/DAt. This overestimates the number of reactions and 
hence decreases the survival probablity, consistent with the results shown in the inset of Fig. [H 
Another way of putting this that in our BD algorithm, the intrinsic association rate is higher than 
that used in the analytical calculations. 



B. Reversible Reactions 



We extend now the dynamical test performed above to the case of the reversible reaction A+B ^ 

C. An analytical solution to the problem, Previf,t\''"Oj^o) is known for one A and one B particle 
[2^. In this test, we adopt the same setup and a similar procedure as for the irreversible case, 
except that we do not stop the run after a reaction, but we let the particle dissociate. At t = tsim we 
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check whether or not the B particle is in the bound state. If it is not, we record the final position. 
The histogram of final positions of the B particles is normalised to the number of survivors at 
^ = ^sim, and compared with the analytical curve Prevl'") ^simko) 0) [2^. The fraction of runs ending 
in the unbound state yields an estimate for the survival probablity 5'rev(isimko)- Again, we initially 
put the B particle at contact (rg = R), so that a larger number reactions and dissociations can 
happen within tgi^. This choice will provide a stringent test for the dynamics of the system. The 
parameters of the system are the same as in Section IIII At with the addition of the dissociation 
rate A;d = lT^^. Similar results were obtained for larger values of k^t^, as well as for other values of 
ro, D, and /ca. 

In Figure [21 we plot prev{r,tsim\R,0) for ^sim ranging from 10~^r to 10~^t, with the same time 
steps as in the previous test. We find the BD algorithm correctly reproduces both the shape and 
the area of the analytical probability distribution. Similarly to Figure [H we show in the Inset Prev, 
computed for tgim = W~^t, for decreasing time steps At: as expected on the basis of the results 
for the irreversible reaction (Fig. [1]), simulations with large time steps underestimate the survival 
probability. 

For the next tests, we consider a setup similar to that used above, but we allow for multiple 
B particles, and we enclose our system in a cubic simulation box of volume V, endowed with 
reflecting walls. All the B particles can bind to the single A particle, and do not interact among 
themselves. This last assumption, valid only in the limit of low packing fractions (p, is satisfied in 
our simulations where (p<0.02. Conversely, particles B and C interact as hard spheres, i.e. they 
are not allowed to overlap. We investigate whether equilibrium properties of the system, such as the 
probability of being in the bound state C (pbound)j are correctly reproduced by the BD algorithm. 
The probability Abound can be evaluated by measuring the time when the C particle is present in 
the system, with respect to the total simulation time. The mean field value for this quantity can 
be obtained from the macroscopic rate equation in steady state: 

Pbound - ^^^^^ ^ , {^V 

where K^q = k^/kd, and V* = V — ^7t{Ra + Rb)^- We note here that this prediction should be 
valid when the radial distribution at contact equals unity; given the low density of particles in our 
system this should be the case. 

We simulate the system with a varying number Nb of B particles, with a fixed time step 
At = W-'^T. We choose K^q = V, so that 

PboundC-^B — 1) — 0.5. The enclosing box measures 
(20i? X 20R X 20R). Figure [3] compares the results of our simulations with the prediction of Eq. [2^ 
we see a clear agreement. To illustrate that obeying the detailed-balance condition is important, we 
also performed a series of simulations in which the particles after dissociation were put at contact; 
in other words, we considered a function g{r) = 5{r — R). This move does violate detailed balance 
and indeed affects a correct estimate for Pbound; as shown in Inset A of Figure [3) the incorrect 
procedure overestimates the time the particle spends in the bound state, especially for a low number 
of B particles. These data clearly show that a naive treatment of the dissociation events leads to 
systematic errors, which are especially severe when the system has a low number of reactants, as 
it is often the case in biochemical networks. Finally, we tested whether the equilibrium properties 
of the system do not depend on the chosen time step. To this end, we compute Pbound 

for Nb = 1 

and different values of the time step At. As illustrated in the Inset B of Figure [3l we obtain a 
good agreement even for very large time steps, where probably the dynamics of the system is not 
entirely natural. 

Finally, we compare our Brownian Dynamics algorithm with the Stochastic Simulation Algo- 
rithm, based on a Kinetic Monte Carlo scheme that propagates the system according to the solution 
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of its zero-dimensional chemical master equation [1^ . This scheme accounts only for the stochastic- 
ity arising from the fluctuations in the number of particles; spatial fluctuations due to the diffusive 
motion of particles are completely neglected. The system is thus assumed to be well-stirred at all 
times. We consider the reversible reaction A+B ^ C for Nb = 1: in the SSA, the association times 
follow a Poisson distribution, with mean l/k{, where = l/{ATTDR) + k^'^ is the macroscopic 
forward rate. 

We collect the association times for a BD run with V = 64000i?^ , At = W^t, = WOR^/t, kos = 
1000r-^ and we compare it with an SSA run obtained with the same set of parameters, but using 
the modified association rate k{. Figure H] compares the two distributions: the BD line shows a 
marked increase in the number of association events at short times, as compared to the Poissonian 
distribution with mean k{ of the SSA. This effect has a purely spatial origin and has been previously 



observed ([181. |22||): when particles dissociate in space, their distance is still very small, therefore 



the probability of an immediate rebinding in next few times steps is very high. Long association 
times, in a BD simulation, are related to particles which have wandered diffusively in the box, and 
have finally found the target. The distribution of such times is again exponential, with a constant 
h. 

The test above show that our BD algorithm, which rigorously obeys detailed balance, correctly 
reproduces the equilibrium properties and provides a good description of the dynamics of the 
system in time and space. 



IV. APPLICATION: THE PUSH-PULL MODEL 

In this Section, we apply our Brownian Dynamics to a simple model of a push-pull network. In 
this network, two antagonistic enzymes continually covalently modify and demodify a substrate, 
respectively; a well-known example is a protein that is phosphorylated and dephosphorylated by 
a kinase and a phosphatase, respectively. The first enzyme converts a substrate molecule into an 
"active" state: bearing in mind the phosphorylation example, we call this active substrate Sp, and 
the enzyme K (kinase) . A molecule in the active state Sp can be brought back to the original state 
S under the action of a second enzyme, P (phosphatase). The model is nicknamed "push-pull", 
as the substrates are continuously switching between the two states, while consuming energy. The 
reactions with the enzymes are described according to Michaelis-Menten kinetics: the two reactants 
form first an intermediate bound state, which can lead either to a dissociation or to the release of 
a converted molecule. In (23 |. the model is solved at the level of the Macroscopical Rate Equation 
at steady state, which yields the average behavior of the system. 

Goldbeter and Koshland showed that such a system can display an ultrasensitive behavior (that 
is, a sensitivity curve steeper than the conventional response showed by the Michaelis-Menten 



mechanism) without the need of introducing cooperative interactions 2j]. More precisely, the 
interplay between two converter enzymes operating in opposite directions on a target whose quan- 
tity is conserved can give rise to a switch-like response in the steady-state fraction of modified 
molecules, when the ratio between the conversion rates is varied. The requirement for such a sharp 
transition is the saturation of the enzymes: the effective conversion rates then become independent 
on the number of substrate molecules, thus making the reaction rates "zero-order" in the substrate 
concentration. 

The above-mentioned analysis does not however account for any kind of fluctuations that may 
arise from the low number of reactants, the stochastic behavior of the chemical reactions, or 
the diffusion of the molecules in space. In [2H], the same model is studied at the level of the 
chemical master equation, taking into account finite-size effects that arise in real systems, that is 
the discreteness and the possible low copy number of enzymes and substrates. In order to achieve 
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ultrasensitivity, the enzymes must be saturated, and therefore their concentration is hkely to be 
very low. Large fluctuations are then observed around their average behavior: the authors show 
that the results obtained with a mesoscopic approach reduce to those of the macroscopic analysis 



of [2j] only when the number of molecules is sufficiently large. If this is not the case, as it can 
easily happen in a bacterial cell where some species are present only in few dozens of copies, the 
sensitivity of the system is reduced, and the response is less steep than the macroscopic theory 
would predict. This deviation can be easily understood when one realises that high sensitivity 
corresponds to highly saturated enzymes. In this regime, the reaction rates do not depend on the 
number of substrate molecules. The system performs a random walk in the number of S molecules 
and it is thus subject to large fluctuations. Our Brownian Dynamics algorithm now allows us to 
study the effect of spatial fluctuations due to the diffusive motion of the molecules. 
The system we consider is defined by the following set of reactions: 

Reaction Rate 

S + K ^ KS fconi, ^offi (25a) 

KS ^K + Sp ki (25b) 

Sp + P^PSp kon2,koS2 (25c) 

PSp^P + S k2. (25d) 

Here k{ stands for the macroscopical association rate. The system will be simulated with the BD 
algorithm in a rectangular box of dimensions Xbox = 20i?, ybox = lOi?, Zbox = lOi?, with a single 
kinase and a single phosphatase molecule, held fixed at distance A = 0.5xbox on the central axis of 
the box, as depicted in Figure [5l The system is initially prepared with A''5tot particles, distributed 
in the two states according to the solution of the macroscopical rate equation. In the following, we 
investigate the effect of spatial fluctuations of the substrate molecules on the input-output relation 
of the system, and we compare the BD results to those obtained with the mean-field and the 
zero-dimensional chemical master equation approach. 

The input-output relation is defined as the mean fraction of phosphorylated substrate molecules 
{Sp) /Nstot as a function of the ratio ki/k2. We compute it with 80 substrate molecules in the 
simulation box, in order to meet the requirement Ns'^Nk and Ns^, Np (Ns + Ns^, = Nstot)- 
The parameters governing the steepness of the sigmoid curves are the Michaelis-Menten equilibrium 
rates Ki = ( fcofn + ) / ^oni and K2 = (A;ofr2 + ^2 ) / ^on2 • In all our simulations we set Ki = K2 = Km ■ 
In our simulations, diffusion constants are in the order of 10~^ — 10~^i?^/At, intrinsic association 
rates between 0.001 and 0.006 R^/At, dissociation rates vary between 5 • 10"'^ and 5 • 10~^(At)~^, 
and production rates are chosen in the range 10"'^ — 10~*(At)~^. To vary ki/k2 keeping Kyi 
constant, we vary /cofn and ki together, keeping their sum constant. When i^M/[5'tot] ^ 1 the 
enzymes are totally saturated and the change in the fraction of modified proteins is abrupt; on 
the other hand, when KM/[Stot\ > 1, the rise of the curve becomes asymptotically close to the 
hyperbolic Michaelis-Menten shape. 

Figure [6] shows several examples of the input-output relation, obtained with three different 
methods: with the analytic macroscopic approach of [i^] (solid lines), with SSA simulations as 
in [2H] (diamonds) and with the Brownian Dynamics algorithm (circles). The BD simulations are 
performed with intrinsic association and dissociation rates k^ and k^, respectively. In the mean- 
field analysis and the SSA simulations of the zero-dimensional chemical master equation, for the 
association reaction the rate constant was chosen to be that of the macroscopic association rate /con, 
as given by fcon = (l/^a+l/fco)"^) with fca being the intrinsic association rate, and fco = '^irRD the 
diffusion- limited rate (see Eq. [5|) . For the rate constant of the backward reaction in the mean-field 
analysis and the SSA simulations, we chose the intrinsic reaction rate and not the macroscopic one 
given by Eq. [6l The reason is that the macroscopic dissociation rate takes into account that, upon 
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dissociation, the dissociated species rebind a number of times before they diffuse away from each 
other into the bulk [3]. As we have recently shown, association and dissociation reactions can be 
described with effective rates given by Eqs. [5]and[6]when the associated species can only dissociate, 
thus when there is no competing decay channel for the associated species [l8(. In particular, in Ref. 
[l^ we studied the effect of spatial fluctuations due to the diffusive motion of repressor molecules 
on the noise in the expression of a gene; the simulation results showed that the stochasticity in 
the binding of the repressor to the DNA resulting from the spatial fluctuations of the repressor 
molecules can be a major source of noise in gene expression; however, this result could be described 
by renormalising the intrinsic association and dissociation rates for repressor-DNA binding using 
the expressions of Eq. [5] and Eq. [6l Here, the situation is markedly different. The reason is 
that the associated species, KS and PSp, can either dissociate or lead to a chemical modification 
reaction, upon which the product must diffuse to the other enzyme in order to be demodified. 
These reaction channels compete with one another, and if they compete with one another on the 
same time scale, the effective dissociation rate is difficult to determine. In our system, however, 
the number of rebindings is in the order of unity. Therefore, renormalising the rate constants does 
not substantially change the actual values of the rate constants. We carried out simulations where 
we chose either to renormalise both the association and dissociation rates or neither of them, and 
we found analogous results. In the following, we show only BD data obtained with the intrinsic 
dissociation and association rate. Ky[ is computed with the intrinsic dissociation rates, and and 

with kon = {k:^^ + K^y^ ■ 

The different panels of Figure [6] show the data for increasing Kyi/lStot]- In panel D, -ftTivi/i'S'tot] = 
3.7 and the response of the system is very similar to a Michaelis-Menten kinetics. In this case, 
the results of the simulations, obtained both with BD and the SSA, perfectly follow the analytical 
curve. When Km > 1, both reactions are in the first-order regime, which means that their rates are 
proportional to the number of substrate molecules. As a result, when this number changes, the rates 
of conversion in the opposite direction changes immediately. This counteracts the modification and 
reduces the effect of fluctuations. However, as we decrease Km, the numerical data start to deviate 
from the predicted macroscopic behaviour. In the case of SSA, the deviation is mild, and barely 
visible in panel A, where i^M/['S'tot] = 0.05. In contrast, the BD simulations yield a much more 
marked deviation, clearly seen already from panel C, where Ku/lStni] = 0.52. 



The data in Figure [6] confirm and extend the findings of Ref. [251 ]: stochastic fluctuations 
of the system dampen the ultrasensitivity, which could be obtained only in an infinitely large, 
well-stirred system. The SSA correctly accounts for the temporal fiuctuations arising from the 
stochastic behavior of chemical reactions, and for the discreteness of the components, but it does 
assume a well-stirred system, where spatial fluctuations can be neglected. These hypothesis result 
in a deviation on the order of few per cents in the ultrasensitive regime. Brownian Dynamics, on 
the contrary, properly accounts for temporal and spatial fluctuations. These last are related to 
the diffusive motion of the substrate molecules; when the concentrations of the species are low, 
they become a serious limiting factor which notably reduces the sensitivity of system. Brownian 
Dynamics is therefore able to show that the response of the system can be much less sensitive than 
predicted at the macroscopic level in the ultrasensitive regime, if the species move slowly enough 
in space, i.e. when the system is not well-stirred. 

Finally, we emphasise that Brownian Dynamics algorithms can be used to directly measure 
spatial properties of the system. Among several possibilities, we choose to show in Figure [7] the 
spatial density of particles along the main axis of the box. Data are obtained for two different values 
of the Michaelis-Menten constant: Ku/iStot] =4.7 (D = IQ-^R^ / /\t,kon = 0.007R^ /At,ki + 
feoffi = 2 • 10-'^{At)-^,Ns,,, = 80), where the system is in the flrst-order regime, and -ftTivi/l'S'tot] = 
0.22 {D = W-^R'^/At,kon = 0.01R^/At,ki + kosi = 5 • 10-^{At)-^,Ns,,, = 80), corresponding 
to the ultrasensitive regime. Substrate molecules are very often bound to the enzymes: at the 
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kinase enzyme, the S density has a sharp peak, and so does the Sp density at the location of the 
phosphatase enzyme. The height of the peak is bigger when the system is simulated for a low value 
of Km, where enzymes are completely saturated. Interestingly, the concentration of both molecules 
shows a gradient along the x direction, higher in the half box where the molecules are produced. 
Such gradient is not particularly appreciable for low values of Km- in this regime, enzymes are 
saturated and dissociation events are rare. The particles have therefore the time to spread in the 
box and reach an homogeneous concentration only occasionally perturbed by a production event. 
For high values of Km , particles are produced more often and they do not have the time to stir in 
the box, leading to an accumulation close to the production sites. The profiles for the two species 
are completely symmetric, as these simulations are obtained for ki = k2- Naturally, such a spatial 
property could not be captured if the system is simulated at the level of the zero-dimensional 
chemical master equation. 



V. DISCUSSION 

In this manuscript, we have presented a Brownian Dynamics algorithm that rigorously obeys 
detailed balance for equilibrium reactions. Consequently, the equilibrium properties of biochem- 
ical networks, such as promoter and receptor occupancies, are reproduced exactly to within the 
statistical error. Moreover, the association and dissociation reaction moves are constructed such 
that they allow for a meaningful interpretation: as the time step At 0, the association and 
dissociation rates approach the intrinsic values corresponding to the reaction rates of the species at 
contact. This is useful, also because it allows the BD results to be compared to theoretical results 
on diffusion-influenced reactions, which describe reactions as reactions between species at contact 



17l |. The results shown in Figs. [T] and [2] show that as At 0, the BD simulations correctly 
describe not on the equilibrium properties, but also the dynamical properties of a bimolecular 
reaction. For larger time steps, the BD results deviate from the analytical results, but this is to 
be expected since the analytical results assume that the molecules move by diffusion up to the 
smallest length and time scales, and that reactions only occur once the molecules have moved by 
diffusion into contact. We believe that while the BD results and the analytical results match for 
At < 10~^i?^/Z) (Figs. [Hand [2]), the BD algorithm gives a good description of the dynamics over 
a large range of time steps, i.e. for At < 10^^, because in this range the BD results can be fitted 
to the analytical results with a different kg, and A;^ (data not shown). 

As an illustrative example, we have applied our BD scheme to a model representing the dy- 
namics of a substrate molecule under the action of two antagonistic enzymes. This model was 



previously analysed with deterministic methods 2j], which revealed an ultrasensitive behavior in 



the response of the system when the enzymes are fully saturated. A study conducted at the level 



of the Chemical Master Equation [25|], thus accounting for the low copy number of the substrate 



molecules, highlighted that the ultrasensitivity predicted in Ref. 2j] cannot be achieved when the 
concentration of the substrate is very low. Temporal fluctuations limit then the sensitivity of the 
system in the ultrasensitive regime. We repeated the analysis of Ref. [2H] simulating the system 
with the SSA, and confirmed their findings. Furthermore, we have investigated the role of spatial 
fluctuations on the system with BD simulations. Our analysis shows that the sensitivity of the re- 
sponse curve in the ultrasensitive regime is furtherly reduced when the diffusion of the reactants is 
taken into account explicitely. In particular, when the diffusion of particles is slow and the system 
is far from well-stirred, spatial fluctuations are the dominant source of noise, and the reduction in 
the gain is significant. 
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APPENDIX A: DETAILED BALANCE FOR A WELL-STIRRED MODEL 

In the case of the well-stirred model we used in Sec. IIIIBI andllV l the detailed-balance condition 
is simpler than in the spatially resolved model. Let Na, Nb, Nc be the number of A, B and C 
molecules and V the volume of the system. The configurational partition function of the system 
can be written as the following sum of terms in the canonical ensemble: 

Q = Y,Q{Na,Nb,Nc), (A1) 

{JV} 

where {A''} denotes all possible combinations of {Na,Nb,Nc}; note that we integrated here over 
the momenta. The choice of the canonical ensemble is motivated by the assumption that the cell 
is a closed system, and it does not exchange particles with the environment. 

Let us consider the case where {A, B, C} are ideal particles in a volume V , except for the fact, 
of course, that A and B can form C. The configurational integral Q for {Na, Nb, Nc} particles is 
then: 

Na Nb Nc 



-,Na ^Nb ^Nc 



'i'A,cm^B,cm5c,cm ^ tNa+ Nb+Nc 



Na\Nb\Nc\ 

where qA is the molecular partition function for an A particle, and the factor 1/{Na^-) takes into 
account the indistinguishability of the A particles. The molecular partition function is given by 
qA = qAQA., ,cm, where QA,cm is the partition function corresponding to the internal degrees of freedom 
relative to the center of mass and qA=V is the partition function associated with the translational 
degrees of freedom of the center of mass. The probability that the system has {Na, Nb, Nc} 
molecules, P{Na, Nb, Nc), is then: 

P{Na, Nb,Nc) = Q{Na, Nb, Nc)/Q. (A3) 

Let us now consider the transition from {Na, Nb, Nc} to {Na — 1, Nb — 1, Nc + 1} molecules. 
The ratio between the probabilities of being in the state after and before the transition is: 

P{Na-1,Nb-1,Nc + 1) NaNb 1 (?c,cm ^^^^ 



P{Na,Nb,Nc) Nc + lV QA^craQB 



cm 



Nc + lV "'^ Nc + IVK 
Please note that -fCgq has dimension of volume, such that the expression on the right-hand-side 



is indeed dimensionless. The above expression serves to illustrate the detailed-balance rule 19|], 
which states that 

-funbound^u^b — -Abound -fb^u- (A6) 

Here Punbound is the probability of being in the state {Na,Nb,Nc}, Pu^h is the probability of a 
transition from {Na, Nb, Nc} to {Na — 1, Nb — 1, Nc + 1}, -Pb-»u is the probability of the reverse 
move, and Pbound is the probability of being in the state {Na — 1, Nb — 1, Nc + 1}. Using Eq. ()A5p 
and the former relation we obtain: 

Pn^h = ^NANB and P^,^^ = K{Nc + I). (A7) 

These transition probabilities precisely correspond to those used in Monte Carlo simulations of the 
zero-dimensional chemical master equation 
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APPENDIX B: DERIVATION OF g{r) 
The function g{r), described in Section |TT] is given by the integral of equation (j20p : 



1 



9{r) 



(vra2)3/2 7o 



r'^dr' / smOdO i dip exp 



„'2 



2rr' cos 6 + r'^ 



(Bl) 



where cj^ = AD At. 

Elementary methods can be used: integration over the angular variables yields 



air) 



1 exp (-r^/fj^) rR 



exp 



r — 2rr' 



exp 



(B2) 



Integrating over all the possible final positions corresponding to an overlap between the particles 
(0<r'<i?), gives 



9{r) 
where 



a_]_ 

vr 2r 



exp 



(r + Rf 



exp 



(r - Rf 



1 




+ 2 


erf ^ 







'r + R\ /-r + R 

+ erf 



(B3) 



erf(x) = 4= Te^'/'dt. 
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LIST OF FIGURES 

1. Radial probability distribution for an irreversible reaction. The four curves refer to different 
tsim and were obtained with time steps At = lO^^tgim) except for tsim = 0.1r where we used 
At = 10~^tsim = 10~^T (r = R^/D, R = Ra + Rb)- Particles were initially positioned at 
contact: rQ = R. The intrinsic association constant is A;a = 1000i?^/r. The numerical results 
(symbols) are in excellent agreement with the analytical curves (solid lines) In the Inset, 
we plot the probability distribution for tsim = 10~^r for several time steps. For large At the 
BD algorithm deviates from the analytical line and underestimates the survival probability. 
Error bars are smaller than symbol sizes. 

2. Radial probability distribution for a reversible reaction. The four curves refer to different 
tsim and were obtained with time steps tstep = 10~^tsimi = R, except for tgim = O.lr 
where we used At = 10~^tsim = 10~^r (r = R^/D, R = Ra + Rb)- Particles were initially 
positioned at contact (ro = i?), and the association rate constant is A;a = 1000iZ^/r {t = R'^/D, 
R = Ra + Rb), while the dissociation rate constant is set to fed = 1t^^. The numerical 
results (circles) agree with the analytical curves (solid lines). In the Inset, the probability 
distribution for tsim = 10~^r is plotted for several values of At: for large values of the time 
step, the BD algorithm deviates from the analytical line and underestimates the survival 
probability. Error bars are smaller than symbol sizes. 

3. Probability of having an A particle bound to a i? particle, as a function of the number of 
B particles. The time step is set to At = 10~^r (r = R? /D, R = Ra + Rb), the intrinsic 
association constant to fcon = 71R^/t so that Pace, f = 0.1. k-^is chosen so that K^q = ki/k\^ = V 
(y = 8000^3) and therefore 

Pbound (-^B — 1) — 0.5. The numerical data obtained with BD are 
in agreement with the mean- field values. The error bars of the numerical results are smaller 
than the size of the circles. In Inset A, the simulations are performed positioning dissociated 
particles at contact. This move violates detailed balance and yield an incorrect Pbound for 
low number of particles. In Inset B, pbound; for ^B = 1 is plotted against the time step 
used in the simulations. To keep Pace f = 0.1, we varied k^n from 22A2R^/t (At = 10-V) to 
0.00026P^/r (At= lO^^r). As expected for an equilibrium quantity, pbound does not depend 
on the chosen time step. 

4. Distribution of association times, for the reaction yl + P <-> C, obtained with the Brownian 
Dynamics algorithm (solid line) and with a Stochastic Simulation Algorithm (dashed line) 
neglecting spatial effects. The data are obtained for y = 64000P^, -D = P^/r, At = 10~^r, fca = 
lOOP^/r, /cd = lOOOr"^ {t = R?/D, R = Ra + Rb)- Spatial simulations account for immediate 
rebindings after a dissociation event, and show a higher probability for short association 
times. The two curves decay exponentially to zero with the 

5. Snapshot of the push-pull system. The pink and the green sphere represent, respectively, 
the kinase and the phosphatase molecule, held fixed along the main axis of the box. The 
blue and red spheres represent S and Sp molecules, respectively. The system is represented 
for Nstot = 50 and a box of 20P x lOP x lOP. 

6. Fraction of converted molecules as a function of ratio of the convertion rates ki/k2- The 
data are shown for increasing values of i^M/['S'tot] from panel A to D. Panel A corresponds 
to full saturation of enzymes, whereas in panel D the system is in the first-order regime. 
The continuous lines are obtained by solving the Macroscopical Rate Equation, whereas 
diamonds correspond to the numerical solutions of the master equation (obtained with the 
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conventional SSA) and circles to the output of our BD simulations. SSA data (error bars are 
smaller than the sizes of the symbols) show a mild deviation only when the system displays 
an ultrasensitive behavior, as in panel A. Brownian Dynamics simulations deviate notably 
from the macroscopic curve when KM/[Stot] < 1- Methods accounting for the stochastic 
behavior of the system show thus a reduction in sensitivity for low values of Kyi- In the 
simulations, we vary ki and fcoff2 keeping their sum constant, so that ki/k2 is varied while 
Km does not change. For all panels, D = 10~^ R'^ /{At), Ns^^^ = 80. Panel A: fci + koSi = 
k2 + koS2 = 5 • 10-6 fca = 0m6R^/{At). Panel B: ki + k^m = k2 + k^m = 5 • lO-^At-^ 
fca = OmiR^/{At), Panel C: ki + kosi = ^2 + fcoff2 = 5 • IQ-^At'^, k^ = 0.003i2^/(At), 
Panel D: ki + k^si = k2 + k^m = 2 • 10"^ Ai'^, k^ = 0.0015i?3/(Ai). 

7. The spatial density profiles for S, Sp {ki/k2 = 1) show clear symmetric peaks around the 
locations of the two enzymes: respectively, the S density is peaked around the kinase enzyme, 
at x = 0.25a;box) and the Sp density around the phosphatase at z = — 0.25a;box- These peaks are 
more pronounced when the system is in the ultrasensitive regime (thinner lines). Moreover, 
for high values of the Michaelis-Menten constants (thicker lines) , the spatial density of the 
particles show a gradients, higher close to the production sites of the molecules. For low 
values of Km the gradients are not appreciable anymore: the particles can diffuse in the box 
and reach an homogeneous distribution, as production events happen on slow time scales. 
Simulation parameters: Ku/lStot] = 4.7,1) = IQ''^ R^ / {At) , N s,,, = 80, fci = k2 = feoffi = 
koS2 = W-^At-\k, = 0.007/?3/(At) Ku/iStot] = 0.22, L» = W-^Ry{At),Ns,,, = 80,ki = 
k2 = feoffi = fcoff2 = 2.5 • W-^At"\ k^ = 0.01R^/{At) 
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